#include<../data_tools/usefull_routines.c>


/*
  This program does the decomposition between disk and bulge in a
  nbody simulations
*/


int main(int argc, char *argv[])
{
  
  //Variables
  int i, j, nbin, type, center;
  int ngas, ndm, ndisk, nsmmp, ntmmp, nfmmp;
  int *indaux, *ind0aux, *ind1aux, *ind2aux, *ind3aux, *ind4aux, *ind5aux;
  int *ind_disk;
  double Jp, Jz, Jdot, Jdotdisk, Ldisk, Eps;
  double v2, E, Ek, Ep;
  char *inputfile, namefiles[1000];
  CM cmbar, cmsat;
  
  /* ./exec <infile> <number_of_epsilon_bins> <flag to activate the
     calculation of the center of the satellite>*/

  inputfile = argv[1];
  nbin      = atoi(argv[2]);
  center    = atoi(argv[3]);
  
  printf("\n");
  printf(" #############################################################\n");
  printf(" Reading Bulge and disk for snapshot:\n");
  printf("%s\n",argv[1]);
  printf(" #############################################################\n");
  printf("\n");
  
  //Load data file
  readGadgetBinary("option_makefile.dat", inputfile);
  
  //Take the particles's id
  indaux = malloc(N_tot*sizeof(int));
  if(indaux == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
  
  for(i=0; i<N_tot; i++) 
    indaux[i] = particles[i].id;
  
  ngas  = Header.npartTotal[0];
  ndm   = Header.npartTotal[1];
  ndisk = Header.npartTotal[2];
  nsmmp = Header.npartTotal[3];
  ntmmp = Header.npartTotal[4];
  nfmmp = Header.npartTotal[5];
  
  ind0aux = malloc(ngas*sizeof(int));
  if(ind0aux == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
  
  ind1aux = malloc(ndm*sizeof(int));
  if(ind1aux == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
  
  ind2aux = malloc(ndisk*sizeof(int));
  if(ind2aux == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
  
  
  
  for(i=(ngas+ndm); i<(ngas+ndm+ndisk); i++) 
    {
      ind2aux[i-ngas-ndm] =  particles[i].id;
    }
  
  
  //Computing center of mass for barions
  printf("Computing center of mass for barions\n");
  centerMass(ind2aux, ndisk, 0, 0, &cmbar);
  compute_angmom(ind2aux, ndisk, 0, 0, &cmbar);
  printf("RCM : %e %e %e\n",cmbar.cm[0],cmbar.cm[1],cmbar.cm[2]);
  printf("VCM : %e %e %e\n",cmbar.vcm[0],cmbar.vcm[1],cmbar.vcm[2]);
  printf("LCM : %e %e %e\n",cmbar.lcm[0],cmbar.lcm[1],cmbar.lcm[2]);
  printf("Done.\n");
  printf("\n");
  
  //Traslating galaxy to center of mass of de barions
  printf("Traslating galaxy to center of mass of de barions\n");
  traslate(&cmbar, indaux, N_tot);
  centerMass(ind2aux, ndisk, 0, 0, &cmdummy);
  compute_angmom(ind2aux, ndisk, 0, 0, &cmdummy);
  
  printf("NEW RCM : %e %e %e\n",cmdummy.cm[0],cmdummy.cm[1],cmdummy.cm[2]);
  printf("NEW VCM : %e %e %e\n",cmdummy.vcm[0],cmdummy.vcm[1],cmdummy.vcm[2]);
  printf("NEW LCM : %e %e %e\n",cmdummy.lcm[0],cmdummy.lcm[1],cmdummy.lcm[2]); 
  printf("Done.\n");
  
  //Orientate the galaxy with the angular momentum of the barions
  printf("\n Orientate galaxy with the angular momentum of the barions\n");
  printf("using %f %f %f\n", cmdummy.lcm[0], cmdummy.lcm[1], cmdummy.lcm[2]);

  float posL[3];

  posL[0] = cmdummy.lcm[0];
  posL[1] = cmdummy.lcm[1];
  posL[2] = cmdummy.lcm[2];
  rotate(posL, indaux, N_tot);
  centerMass(ind2aux, ndisk, 0, 0, &cmdummy);
  compute_angmom(ind2aux, ndisk, 0, 0, &cmdummy);
  
  printf("NEW RCM : %e %e %e\n",cmdummy.cm[0],cmdummy.cm[1],cmdummy.cm[2]);
  printf("NEW VCM : %e %e %e\n",cmdummy.vcm[0],cmdummy.vcm[1],cmdummy.vcm[2]);
  printf("NEW LCM : %e %e %e\n",cmdummy.lcm[0],cmdummy.lcm[1],cmdummy.lcm[2]); 
  printf("Done.\n");
  printf("\n");
  
  //Printing galaxy in the certer of mass of the barions
  printf("Printing galaxy in the certer of mass of the barions in *.cmb.*\n");
  sprintf(namefiles,"%s.cmb",inputfile);
  printAscii(namefiles,0,-1);
  printf("Done.\n");
  printf("\n");
  
  
  //changing satellites's particles in its center of mass
  ind1aux = malloc(ndm*sizeof(int));
  if(ind1aux == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  for(i=ngas; i<(ngas+ndm); i++)
    { 
      ind1aux[i-ngas] =  particles[i].id;
    }
  centerMass(ind1aux, ndm, 0, 0, &cmsat);
  compute_angmom(ind1aux, ndm, 0, 0, &cmsat);
  traslate(&cmsat, ind1aux, ndm);
  free(ind1aux);
  
  ////////////////////////////////////////////////////
  /*          This is done for the satellites       */
  
  ind3aux = malloc(nsmmp*sizeof(int));
  if(ind3aux == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  for(i=(ngas+ndm+ndisk); i<(ngas+ndm+ndisk+nsmmp); i++)
    { 
      ind3aux[i-ngas-ndm-ndisk] =  particles[i].id;
    }
  
  ind4aux = malloc(ntmmp*sizeof(int));
  if(ind4aux == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
  
  for(i=(ngas+ndm+ndisk+nsmmp); i<(ngas+ndm+ndisk+nsmmp+ntmmp); i++)
    { 
      ind4aux[i-ngas-ndm-ndisk-nsmmp] =  particles[i].id;
    }
  
  ind5aux = malloc(nfmmp*sizeof(int));
  if(ind5aux == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
  
  for(i=ngas+ndm+ndisk+nsmmp+ntmmp; i<(ngas+ndm+ndisk+nsmmp+ntmmp+nfmmp); i++) 
    {
      ind5aux[i-ngas-ndm-ndisk-nsmmp-ntmmp] =  particles[i].id;
    }
  
  
  //Printing galaxies in its own center of mass
  printf("Printing galaxies in its own center in *.cmb.ind*\n");
  sprintf(namefiles,"%s.cmb.ind",inputfile);
  for(i=3; i<6; i++)
    {
      if(Header.npartTotal[i]!=0)
	{
	  printAscii(namefiles,0,i);
	}
    }
  printf("Done.\n");
  printf("\n");
  
  
  //Selecting SMMP with bound particles
  if(Header.npartTotal[3]!=0)
    {
      type = 3;
      if(center==0)
	{
	  Energybinding(type, &cmsat);
	  traslate(&cmsat, ind3aux, nsmmp);
	}
      else
	{
	  int *ind_bound;
	  float a,b,c,d;
	  ind_bound = (int *)malloc(100*sizeof(int));
	  
	  FILE *bound_part;
	  bound_part = fopen("bound_particles.3","r");
	  if(bound_part==NULL) printf("No existe el archivo bound_particles.3\n");
	  
	  for(i=0; i<100; i++)
	    {
	      fscanf(bound_part,"%d %f %f %f %f\n",&ind_bound[i],&a,&b,&c,&d);
	    }
	  fclose(bound_part);
	  
	  centerMass(ind_bound, 100, 0, 0, &cmsat);
	  traslate(&cmsat, ind3aux, nsmmp);
	  
	  free(ind_bound);
	}  
      
      printf("\nparticles type %d before selection %d\n",type,nsmmp);
      printf("particles disk before selection %d\n",ndisk);
      
      for(i=0; i<N_tot; i++)
	{
	  if(particles[i].type==type)
	    {
	      
	      v2 = particles[i].vel[0]*particles[i].vel[0]
		+particles[i].vel[1]*particles[i].vel[1]
		+particles[i].vel[2]*particles[i].vel[2];
	      
	      Ek = 0.5*particles[i].mass*v2;
	      Ep = particles[i].pot*particles[i].mass;  
	      E = Ek + Ep;
	      
	      if(E >= 0.0)
		{
		  particles[i].type = 2;
		  nsmmp = nsmmp - 1;
		  ndisk = ndisk + 1;
		}
	    }
	}
      
      printf("particles type %d after selection %d\n",type,nsmmp);
      printf("particles disk after selection %d\n\n",ndisk);
    }
  
  //Selecting TMMP with bound particles
  if(Header.npartTotal[4]!=0)
    {
      type = 4;
      if(center==0)
	{
	  Energybinding(type, &cmsat);
	  traslate(&cmsat, ind4aux, ntmmp);
	}
      else
	{
	  int *ind_bound;
	  float a,b,c,d;
	  ind_bound = (int *)malloc(100*sizeof(int));
	  
	  FILE *bound_part;
	  bound_part = fopen("bound_particles.4","r");
	  if(bound_part==NULL) printf("No existe el archivo bound_particles.3\n");
	  
	  for(i=0; i<100; i++)
	    {
	      fscanf(bound_part,"%d %f %f %f %f\n",&ind_bound[i],&a,&b,&c,&d);
	    }
	  fclose(bound_part);
	  
	  centerMass(ind_bound, 100, 0, 0, &cmsat);
	  traslate(&cmsat, ind4aux, nsmmp);
	  
	  free(ind_bound);
	}  
      
      printf("\nparticles type %d before selection %d\n",type,ntmmp);
      printf("particles disk before selection %d\n",ndisk);
      
      for(i=0; i<N_tot; i++)
	{
	  if(particles[i].type==type)
	    {
	      
	      v2 = particles[i].vel[0]*particles[i].vel[0]
		+particles[i].vel[1]*particles[i].vel[1]
		+particles[i].vel[2]*particles[i].vel[2];
	      
	      Ek = 0.5*particles[i].mass*v2;
	      Ep = particles[i].pot*particles[i].mass;  
	      E = Ek + Ep;
	      if(E>=0.0)
		{
		  particles[i].type=2;
		  ntmmp = ntmmp -1;
		  ndisk = ndisk +1;
		}
	    }
	}
      
      printf("particles type %d after selection %d\n",type,ntmmp);
      printf("particles disk after selection %d\n\n",ndisk);
    }
  
  //Selecting FMMP with bound particles
 if(Header.npartTotal[5]!=0)
   {
     type = 5;
     if(center==0)
       {
	 Energybinding(type, &cmsat);
	 traslate(&cmsat, ind5aux, nfmmp);
       }
     else
       {
	 int *ind_bound;
	 float a,b,c,d;
	 ind_bound = (int *)malloc(100*sizeof(int));
	 
	 FILE *bound_part;
	 bound_part = fopen("bound_particles.5","r");
	 if(bound_part==NULL) printf("No existe el archivo bound_particles.3\n");
	 
	 for(i=0; i<100; i++)
	   {
	     fscanf(bound_part,"%d %f %f %f %f\n",&ind_bound[i],&a,&b,&c,&d);
	   }
	 fclose(bound_part);
	 
	 centerMass(ind_bound, 100, 0, 0, &cmsat);
	 traslate(&cmsat, ind5aux, nfmmp);
	 
	 free(ind_bound);
       }  
     
     printf("\nparticles type %d before selection %d\n",type,nfmmp);
     printf("particles disk before selection %d\n",ndisk);
     
     for(i=0; i<N_tot; i++)
       {
	 if(particles[i].type==type)
	   {
	     
	     v2 = particles[i].vel[0]*particles[i].vel[0]
	       +particles[i].vel[1]*particles[i].vel[1]
	       +particles[i].vel[2]*particles[i].vel[2];
	     
	     Ek = 0.5*particles[i].mass*v2;
	     Ep = particles[i].pot*particles[i].mass;  
	     E = Ek + Ep;
	     if(E>=0.0)
	       {
		 particles[i].type=2;
		 nfmmp = nfmmp -1;
		 ndisk = ndisk +1;
	       }
	   }
       }
     
     printf("particles type %d after selection %d\n",type,ntmmp);
     printf("particles disk after selection %d\n\n",ndisk);
   }

 
 //Computing properties of the particles
 printf("Computing properties of the particles\n");
 particlesProperties(N_tot, 1, center);  
 printf("Done.\n");
 printf("\n");
 
 
 sprintf(namefiles,"%s.cmb",inputfile);
 
 //todo vuelve al cm de los bariones

 //histogram for epsilon with all particles of disk
 
 epsilonHistogram(inputfile,-1,2,nbin);
 
 
 //Selecting type of particles
 
 double angle,Jmax;
 angle = atof(argv[4]);

  
  //Selecting Thin disk with Eps*Je/Jp>=cos(inclination) 
  for(i=0; i<N_tot; i++)
    {
      if(particles[i].type==2)
	{
	  Eps = particles[i].epsilon;
	  Jp = sqrt(pow(particles[i].mang[0],2)
		     +pow(particles[i].mang[1],2)
		     +pow(particles[i].mang[2],2));
	  Jmax = sqrt(pow(particles[i].je[0],2)
		       +pow(particles[i].je[1],2)
		       +pow(particles[i].je[2],2));
	 
	  Jdot =  particles[i].je[0]*particles[i].mang[0]
	    +particles[i].je[1]*particles[i].mang[1]
	    +particles[i].je[2]*particles[i].mang[2];
	 
	  Jz = particles[i].mang[2];

	  Jdotdisk =  cmbar.lcm[0]*particles[i].mang[0]+
	    cmbar.lcm[1]*particles[i].mang[1]+
	    cmbar.lcm[2]*particles[i].mang[2]; 
	  
	  Ldisk = sqrt(cmbar.lcm[0]*cmbar.lcm[0]+
		       cmbar.lcm[1]*cmbar.lcm[1]+
		       cmbar.lcm[2]*cmbar.lcm[2]);
	  
	  if((particles[i].r <= 10.0))
	    {
	      
	      //if(((Eps*Jdot)/(Jz*Jp))>=0.605)  
	      //if(((Eps*Jdot)/(Jz*Jp))>=0.573 && ((Jdotdisk/(Ldisk*Jp))>=0.573)) //mejor para 0.2
	      //if(((Eps*Jdot)/(Jz*Jp))>=angle && ((Jdotdisk/(Ldisk*Jp))>=angle)) //mejor para 0.2
	      //if((Jdot/(Jmax*Jp))>=angle && ((Jdotdisk/(Ldisk*Jp))>=angle))  // ALL
	      //if((Jdotdisk/(Ldisk*Jp))>=angle) // disk
	      //if(((Eps*Jdot)/(Jz*Jp))>=0.786 && ((Jdotdisk/(Ldisk*Jp))>=0.786)) //mejor para 0.7
	      //if((Jdotdisk/(Ldisk*Jp))>=angle)  

	      if((Jdot/(Jmax*Jp)) >= angle) //orbit
		{
		  particles[i].type = 2;
		}
	      else
		{
		  particles[i].type = 6;
		}
	    }
	  else
	    {
	      particles[i].type = 8; 
	    }
	}
    }



  //************************
  //Histogram for thin disk
  //************************
  epsilonHistogram(inputfile, 1, 2, nbin);

  //********************************
  //Histogram for bulge + thick disk
  //********************************

  epsilonHistogram(inputfile, 1, 6, nbin);

  //==================================
  //Rotation curves
  //==================================
  printf("===============\n");
  printf("Rotation curves\n");
  printf("===============\n");



  int N2 = 0;

  for(i=0;i<N_tot;i++)
    {
      if(particles[i].type==2)
	{
	  N2 = N2 + 1;
	}
    }

  printf("\n***disco final con %d particulas***\n",N2);      


  ind_disk = malloc(N2*sizeof(int));
  if(ind_disk == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  j = 0;
  for(i=0; i<N_tot; i++)
   {
     if(particles[i].type==2)
       {
	 ind_disk[j] =  particles[i].id;
	 j = j+1;
       } 
   }
  
  //Computing center of mass for disk
  printf("Computing center of mass for disk\n");
  centerMass(ind_disk, N2, 0, 0, &cmbar);
  compute_angmom(ind_disk, N2, 0, 0, &cmbar);
  printf("RCM : %e %e %e\n",cmbar.cm[0],cmbar.cm[1],cmbar.cm[2]);
  printf("VCM : %e %e %e\n",cmbar.vcm[0],cmbar.vcm[1],cmbar.vcm[2]);
  printf("LCM : %e %e %e\n",cmbar.lcm[0],cmbar.lcm[1],cmbar.lcm[2]);
  printf("Done.\n");
  printf("\n");
  
  //Traslating galaxy to center of mass of de barions
  printf("Traslating galaxy to center of mass of de barions\n");
  traslate(&cmbar, indaux, N_tot);
  centerMass(ind_disk, N2, 0, 0, &cmdummy);
  compute_angmom(ind_disk, N2, 0, 0, &cmdummy);
  
  printf("NEW RCM : %e %e %e\n",cmdummy.cm[0],cmdummy.cm[1],cmdummy.cm[2]);
  printf("NEW VCM : %e %e %e\n",cmdummy.vcm[0],cmdummy.vcm[1],cmdummy.vcm[2]);
  printf("NEW LCM : %e %e %e\n",cmdummy.lcm[0],cmdummy.lcm[1],cmdummy.lcm[2]); 
  printf("Done.\n");
  
  //Orientate the galaxy with the angular momentum of the barions
  printf("\n Orientate galaxy with the angular momentum of the barions\n");
  printf("using %f %f %f\n", cmdummy.lcm[0], cmdummy.lcm[1], cmdummy.lcm[2]);
  posL[0] = cmdummy.lcm[0];
  posL[1] = cmdummy.lcm[1];
  posL[2] = cmdummy.lcm[2];
  rotate(posL, indaux, N_tot);
  centerMass(ind_disk, N2, 0, 0, &cmdummy);
  compute_angmom(ind_disk, N2, 0, 0, &cmdummy);
  
  printf("NEW RCM : %e %e %e\n",cmdummy.cm[0],cmdummy.cm[1],cmdummy.cm[2]);
  printf("NEW VCM : %e %e %e\n",cmdummy.vcm[0],cmdummy.vcm[1],cmdummy.vcm[2]);
  printf("NEW LCM : %e %e %e\n",cmdummy.lcm[0],cmdummy.lcm[1],cmdummy.lcm[2]); 
  printf("Done.\n");
  printf("\n");

  
  
  //Computing properties of the particles
  printf("Computing properties of the particles\n");
  particlesProperties(N_tot,0,center);  
  printf("Done.\n");
  printf("\n");


  //===============================
  //    INCLINATIONS DISK
  //===============================

  //inclination();


  //count particles by type 
  int N[9]={0,0,0,0,0,0,0,0,0};

  ngas  = Header.npartTotal[0];
  ndm   = Header.npartTotal[1];
  ndisk = Header.npartTotal[2];
  nsmmp = Header.npartTotal[3];
  ntmmp = Header.npartTotal[4];
  nfmmp = Header.npartTotal[5];
 
  for(i=0;i<N_tot;i++)
    {
      if(particles[i].type==0)
	{
	  N[0] = N[0] + 1;
	}      
      if(particles[i].type==1)
	{
	  N[1] = N[1] + 1;
	}
      if(particles[i].type==2)
	{
	  N[2] = N[2] + 1;
	}
      if(particles[i].type==3)
	{
	  N[3] = N[3] + 1;
	}
      if(particles[i].type==4)
	{
	  N[4] = N[4] + 1;
	}
      if(particles[i].type==5)
	{
	  N[5] = N[5] + 1;
	}
      if(particles[i].type==6)
	{
	  N[6] = N[6] + 1;
	}
      if(particles[i].type==7)
	{
	  N[7] = N[7] + 1;
	}
      if(particles[i].type==8)
	{
	  N[8] = N[8] + 1;
	}
  
    }

  //sprintf(namefiles,"%s",inputfile);
  //rotation_curve(1,N[1],namefiles);
  //rotation_curve(2,N[2],namefiles);
  //rotation_curve(6,N[6],namefiles);
 // total_rotation_curve(namefiles);
  
 
  sprintf(namefiles,"%s.halo",inputfile);
  printType(namefiles,1);

  sprintf(namefiles,"%s.disk",inputfile);
  printType(namefiles,2);

  sprintf(namefiles,"%s.bulge",inputfile);
  printType(namefiles,6);

  sprintf(namefiles,"%s",inputfile);
  printAscii(namefiles,1,-1);

 
  
  printf("\n*** Angle = %lf ***\n",angle);

  printf("\n====particulas iniciales====\n");
  printf("tipo 1: %d tipo 2: %d tipo 3: %d: Total: %d\n\n",
	 Header.Npart[1],
	 Header.Npart[2],
	 Header.Npart[3],
	 Header.Npart[1]+Header.Npart[2]+Header.Npart[3]);

  printf("\n====particulas finales====\n");
  printf("tipo 0: %d\n",N[0]);
  printf("tipo 1: %d\n",N[1]);
  printf("tipo 2: %d\n",N[2]);
  printf("tipo 3: %d\n",N[3]);
  printf("tipo 4: %d\n",N[4]);
  printf("tipo 5: %d\n",N[5]);
  printf("tipo 6: %d\n",N[6]);
  printf("tipo 7: %d\n",N[7]);
  printf("tipo 8: %d\n",N[8]);
  printf("Total : %d\n",N[0]+N[1]+N[2]+N[3]+N[4]+N[5]+N[6]+N[7]+N[8]);
 
  free(ind_disk);

  return 0;
}
